
__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}

__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	int o = 0;
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
		for (int j=0; j<=ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<=nz; k++)
				f1[i][j] += func[o++] * nn[k];
		}

		GetBernstein(nn, ny, v);
		for (int i=0; i<=nx; i++)
		{
			f2[i] = 0.0f;
			for (int j=0; j<=ny; j++)
				f2[i] += f1[i][j]*nn[j];
		}

		GetBernstein(nn, nx, u);
		(*re) = 0.0f;
		for (int i=0; i<= nx; i++)
		{
			(*re) += f2[i]*nn[i];
		}
}
